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Abstract 

Chaos and oscillations continue to capture the interest of both the scientific and pub- 
lic domains. Yet despite the importance of these qualitative features, most attempts at 
constructing mathematical models of such phenomena have taken an indirect, quantitative 
approach, e.g. by fitting models to a finite number of data-points. Here we develop a 
qualitative inference framework that allows us to both reverse engineer and design systems 
exhibiting these and other dynamical behaviours by directly specifying the desired charac- 
teristics of the underlying dynamical attractor. This change in perspective from quantitative 
to qualitative dynamics, provides fundamental and new insights into the properties of dy- 
namical systems. 

Mathematical modelling requires a combination of experimentation, domain knowledge 
and, at times, a measure of luck. Beyond the intrinsic challenges of describing complex and 
complicated phenomena, the difficulty resides at a very fundamental level with the diversity 
of models that could explain a given set of observations. This is a manifestation of the so- 
called inverse problem (TJ, which is encountered whenever we aim to reconstruct a model of 
the process from which data have been generated. Exploring the potential space of solutions 
computationally can be prohibitively expensive and will generally require sophisticated nu- 
merical approaches or search heuristics, as well as expert guidance and manual interventions. 
Parameter estimation [2], model inference |3j and model selection (4|[5j all address aspects of 
this problem. 

The inverse problem also applies in a different context: the design of systems with specified 
or desired outputs. Here again we have a multitude of different models — or for sufficiently 

*The authors would like to dedicate this article to the memory of Jaroslav Stark. 
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Figure 1: Encoding and inferring the desired dynamics, (a) Lyapunov exponents (LEs), Ao, ...A n , characterise 
the contraction/expansion of an initially small perturbation, eo, to the system, (b) The leading LE determines the 
principal dynamics and characteristics of the attractor of a dynamical system. For Ao < the attractor will be a 
stable fixed-point; stable oscillating solutions will be obtained if Ao = 0; for Ao > we observe chaos and the 
system will exhibit a so-called strange attractor; if more than one LE is positive, then we speak of hyperchaos and 
the attractor will exhibit behaviour with similar statistical properties to white noise, (c) Key steps in the unscented 
Kalman filter (UKF) for qualitative inference. At the k th iteration, the current prior parameter distribution is 
formed by perturbing the previous posterior, Ok, with the process noise Vk- The distribution of LEs for the model 
/ induced by the prior parameter distribution is calculated via the LE estimation routine L and the unscented 
transform. Comparing the mean LE, A&, to the target LE, \ ta rget, the prior parameters are updated using the UKF 
update equations. As the filter proceeds, parameters are found that locally minimise the sum of squared error 
between target and estimated LEs. 

complicated models a potentially vast range of parameters — that fulfil a given set of design 
objectives. Therefore system design can be fraught with the same challenges as statistical in- 
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ference or reverse engineering tasks: in the former case we want to learn the existing structure 
and properties of a system that has produced certain types of data, while in the latter we want to 
design constructible systems that will reliably and robustly exhibit certain types of behaviour. 

These challenges are often further exacerbated by unsuitable or insufficient encoding of the 
behaviour that we observe (in natural systems) or would like to see (in designed systems). For 
example, if we aim to estimate parameters describing an oscillating system from a series of 
observations, then it is possible to get good and even globally optimal fits to the data, with- 
out finding a qualitatively acceptable solution. Various methods of qualitative inference have 
been developed to address this issue; the topology of bifurcation diagrams |6[[7j, local stability 
properties of dynamically invariant sets (8}|9j ? symbolic sequences of chaotic systems JT0| | and 



temporal logic constraints |TT}[T2J have variously been used to drive parameter searches, or for 
model checking. However these methods are either limited in the complexity of behaviour they 
can detect, or by conditioning on surrogate data (e.g. forcing solutions through a small number 
of points), they suffer in the same way as quantitative approaches. The method proposed here 
extends the scope of the promising but underdeveloped class of qualitative parameter estima- 



tion algorithms JT3J, allowing detection and control of the most complex and elusive dynamical 
behaviours, such as oscillations, chaos and hyperchaos. 
We consider models of the general form 

d ^ = f(y(t),y ;e), (1) 

where y(t) denotes the n-dimensional state of the system at time t, f is the gradient field char- 
acterised by a parameter vector, 0, and y = y(0) are the initial conditions, which may be 
unknown, too. Coaxing the solutions of such systems into exhibiting a desired dynamical be- 
haviour is reliant upon the ability to, firstly, encode the behaviour sufficiently as constraints 
upon a set of model properties that may be conveniently evaluated, and secondly, to identify 
regions in parameter space for which these constraints are satisfied. Here we meet these chal- 
lenges using a combination of statistical and dynamical systems techniques. In particular, we 
pose the problem within a state-space framework, where the observation function corresponds 
to evaluating the type of attractor exhibited by the model with given parameters and initial con- 
ditions. We then exploit the flexibility and efficiency of the unscented Kalman filter (UKF) to 
systematically move in parameter space until the desired or expected dynamical behaviour is 
exhibited. The approach, outlined in Fig. 1 and developed fully in the Methods and online 
Supplementary Information, is demonstrated below within different contexts, covering some 
classical dynamical model systems and electronic circuits that exhibit oscillations, chaos and 
hyperchaos, and a biological regulatory system that exhibits oscillatory behaviour. 



Results 

Oscillations and chaos in electronic circuits 

The elimination of chaos from a system, or conversely its "chaotification", have potential ap- 
plications to biological, medical, information processing and other technological systems p4| . 
Here, we use a simple electric circuit [15] (shown in Fig. 2a), to illustrate how our method 
can be used to tune the system parameters such that the dynamics are driven into and out of 
chaos. The circuit model includes a parameter a, representing the scaled resistance of a variable 
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Figure 2: Detecting and controlling chaos and oscillations. Plots show the estimated parameters at successive 
iterations of the unscented Kalman filter. Snapshots of the developing attractor are shown above the plots. The 
sum of squared error (E) is indicated for different sections of the parameter trajectories, (a, b) The filter is able 
to drive the electric circuit between oscillations and chaos in less than 10 iterations, (c) Parameter trajectories 
for a simple model of the Hesl regulatory system that yield oscillations. Several regions in parameter space can 
be identified that exhibit oscillatory behaviour, (d) Examples of trajectories generated from a region in parameter 
space which was found using our approach. Here we used the qualitative inference procedure in order to elicit 
a prior to be used for parameter inference. Trajectories4vere sampled from the prior. Data are indicated by red 
circles and represent fold change in Hesl mRNA; the blue strips indicate the confidence intervals obtained using 
Gaussian Process regression, in which standard Gaussian noise is assumed, with maximum marginal likelihood 
estimates for the other hyperparameters pl| . 



resistor, i? 2 , which we make the lone subject of the inference. In turn we start the system in 
an oscillatory regime and tune the parameter according to the posterior predictions at each step 
of the UKF, until we enter a chaotic regime, and vice versa (see Fig. 2b). The two desired 
behaviours are encoded as constraints only upon the target maximal Lyapunov exponent (LE), 
specifying, Ai = 0, for oscillations, and, Ai = d > S to i for chaos, where S to i is taken larger than 
the expected error in the LE estimation procedure, as discussed in the Supplementary Informa- 
tion online. 

For systems of this size, the qualitative dynamical regimes can be explored exhaustively and 
in short time (finding the desired behaviour takes minutes even for moderate sized systems). 

Detecting oscillations in immune signalling 

Oscillations appear to be ubiquitous in nature, yet for reasons noted above they often remain 
elusive to quantitatively driven parameter inference techniques. Here we consider a dynamical 
system describing the expression levels of the transcription factor Hesl, which is involved in 
regulating the segmentation of vertebrate embryos (16). Oscillations of Hesl expression levels 
have been observed in vitro in mouse cell lines, and reproduced using various modelling ap- 
proaches, including continuous deterministic delay (16} [T7J and discrete stochastic delay mod- 
els JT8) . We investigate a simple three component ODE model of the regulatory dynamics with 
mRNA transcription modelled by a Hill function, 



M = -k deg M + 1/(1 + (P 2 /P ) h ) (2) 
A = -k deg P x + vM-k x P x (3) 
P2 = -k deg P2 + hP 1: (4) 

where state variables M, P 1 and P 2 , are the molecular concentrations of Hesl mRNA, cyto- 
plasmic and nuclear proteins respectively. The parameter k deg is the Hesl protein degradation 
rate which we assume to be the same for both cytoplasmic and nuclear proteins, k\ is the rate 
of transport of Hesl protein into the nucleus, P is the amount of Hesl protein in the nucleus 
when the rate of transcription of Hesl mRNA is at half its maximal value, v is the rate of trans- 
lation of Hesl mRNA, and h is the Hill coefficient. For the inference we take, k deg , to be the 
experimentally determined value of 0.03 mm -1 Jl9| ]. 

In Fig. 2c we show the results for the inference using our algorithm on the model shown 
above. Note that the value inferred for parameter is significantly lower than the range of 
values investigated for the continuous deterministic delay model of H. Momiji and N. A. M. 



Monk p7| |. Interestingly, repeating the inference with different initial parameter sets leads to 
similar values of k\ (k\ < 0.01), but to a broad range of values for the other parameters, all of 
which result in oscillatory behaviour. Our qualitative inference thus suggests that oscillations of 
Hesl protein and mRNA levels are strongly dependent upon maintaining a low rate of transport 
of Hesl protein into the nucleus, and that the dependence on other system parameters is less 
strong. As 1/ki is the expected time Hesl spends in the cytoplasm, this corresponds to the 
delay that had previously been posited to be necessary for such oscillations to occur [17]. Our 
approach readily identifies a parameter regime exhibiting oscillatory dynamics without explic- 
itly requiring (discrete) time-delays. 

Next, we used the qualitative inference result as the basis to estimate the model parameters 
from the Hesl data described below. An approximate Bayesian computation algorithm (ABC 
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Figure 3: Designing attractive models, (a) Inferring a complete spectrum. After only 22 iterations, the charac- 
teristic "butterfly" strange attractor emerges. The final parameters and LEs are a = 10.2, p = 29.2, f3 — 2.45 
and (0.899, 2.74 x 10 -4 , —14.6). (b) A function of the Lyapunov exponents, the Kaplan- Yorke fractal dimension 
may also be used to specify the desired attractor. Here parameters for a target dimension of 1 are found for the 
Lorenz system within 20 iterations, giving rise to a limit cycle as required by the theory, (c) 3 -dimensional projec- 
tions of the hyperchaotic system with parameter vector (a, 6, c, <i, e, /)=(49.98, 35.86, 30.5, 1.35, 36.6, 33.8) and 
corresponding LEs (31.8, 16.8, —19.1,— 71.4). A very chaotic attractor. Within few iterations our algorithm was 
able to drive the system towards an attractor characterised by Lyapunov exponents twice the size of any that had 
previously been reported. 

SMC (20)), capable of sampling from non-Gaussian and multimodal posteriors, was employed 
and Fig. 2d shows the fits of simulated trajectories for 20 parameters drawn randomly from 
the resulting posterior distribution; these are in good agreement with the confidence intervals 
(the blue bands in Fig. 2d), which can be obtained from the time-course data via a Bayesian 
nonparametric method pl| . It is worth noting that using the UKF alone, we could in principle 
consider the Lyapunov exponents and data together in order to infer parameters that are both 
qualitatively and quantitatively acceptable. However by splitting the inference, we take advan- 
tage of the strengths of each algorithm within the Bayesian framework; first we exploit the 
efficiency of the unscented Kalman filter to work with a sophisticated encoding of the desired 
behaviour that is computationally expensive to calculate; subsequently we use this qualitative 
information in order to construct suitable priors for an ABC method capable of dealing with 
non-Gaussian posteriors. 
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Designing attractors 



While the maximal LE alone is sufficient to encode fixed points, limit cycles and strange attrac- 
tors, we may include additional target exponents to design the complete Lyapunov spectrum 
(Fig. 3a), design the (Kaplan- Yorke) fractal dimension |22| (D = k + Y^l=i ^i/|^fc+i| ? where k 
is the largest integer for which ^2 i=1 Xi > 0) of a system's attractor (Fig. 3b), or drive models 
to behave hyper-chaotically (Fig. 3c). 

The first two of these applications are illustrated with the Lorenz system (23) which has 
become a canonical example of how sensitivity to initial conditions can give rise to unpre- 
dictable behaviour. The model is known to exhibit a chaotic regime with Lyapunov exponents, 
A = (0.906,0,-14.57), for parameter vector (a, = (10,28,8/3). Here we infer back 
these parameters, starting with different prior means, by setting our target Lyapunov spectrum 
to A. If we restrict the parameter search to the region [0, 30] 3 , as described in the Supplemen- 
tary Information, we are able to do this reliably from random starting positions. The parameter 
trajectories and evolving attractor of a representative run of the inference algorithm is shown in 
Fig. 3a, where the sum of squared error between estimated and target LEs is less than 8 x 10" 5 
after the 100th iteration. Without constraints on the parameters the inference algorithm con- 
verges to different parameter combinations that display indistinguishable LEs. This allows us 
to assess (for example) the robustness of chaotic dynamics by mapping systematically the re- 
gions of parameter space that yield similar Lyapunov exponents. Fig. 3b shows how the fractal 
dimension - a function of the Lyapunov exponents - may also be tuned (in this example, halved). 
While computational difficulties have in the past precluded such investigations, our approach al- 
lows us to map attractor structures (and the range of parameters giving rise to similar attractors) 
very efficiently. 

For the third application of driving a system into hyper-chaos, we investigate a four di- 
mensional system with six parameters, whose significance lies in having two very large LEs 
(Ai e [10.7741,12.9798] and A 2 e [0.4145, 2.6669]) over a broad parameter range (24). The 
resulting highly complex deterministic dynamics share statistical properties with white noise, 
making it attractive for engineering applications such as communication encryption and random 
number generation. By setting large target values of Ai and A 2 , we use our method to obtain pa- 
rameters for which the system displays LEs that are over two times bigger than previously found 
for the system. Fig. 3c, shows the three dimensional projections of the resulting hyper-chaotic 
attractor. 

These are, of course, toy-applications, but they demonstrate the flexibility and potential uses 
of this approach: we can really start to explore qualitative behaviour in a numerically efficient 
and speedy manner. For example, it becomes possible to map (or design) the qualitative charac- 
teristics of complex systems and to test robustness of qualitative system features systematically. 



Discussion 

We have demonstrated that it is possible to use statistical inference techniques in order to con- 
dition dynamical systems on observed (biological oscillations in Hesl) or desired qualitative 
characteristics (oscillations, chaos and hyper-chaos in natural and engineered systems). This 
provides us with unprecedented ability to probe the workings of dynamical systems. Here we 
have only used the approach for inference and design of the attractors of dynamical systems, as 
encoded by their Lyapunov exponents. This, however, has already been enough to show that it 
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is not necessary to impose discrete time-delays in order to explain the oscillations in the Hesl 
system (T^{T5). 

A focus on qualitative features has several advantages: first, it is notoriously difficult, for ex- 
ample, to ensure that parameter inference preferentially (let alone exclusively) explores regions 
in parameter space that correspond to the correct qualitative behaviour such as oscillations. This 
is the case for optimisation as well as the more sophisticated estimation procedures. Arguably, 
however, solutions which display the correct qualitative behaviour are more interesting than 
those which locally minimise some cost-function in light of some limited data. Obviously, in a 
design setting ensuring the correct qualitative behaviour is equally important. 

Second, the numerical performance of the current approach allows us to study fundamental 
aspects related to the robustness of qualitative behaviour. This allows us for the first time to 
ascertain how likely a system it is to produce a given Lyapunov spectrum (and hence attractor 
dimension) for different parameter values, 9. Our approach, coupled with means of covering 
large-dimensional parameter spaces, such as Latin-hypercube or Sobol sampling J25] l, allows 
us to explore such qualitative robustness. Or more specifically, we can map out boundaries 
separating areas in phase space with different qualitative types of behaviour. We can also drive 
systems into regions with Lyapunov exponents of magnitudes not previously observed. The last 
aspect will have particular appeal to information and communication scientists as such hyper- 
chaos shares important properties with white noise and potential applications in cryptography 



and coding theory abound [26]. 



Finally, our approach can also be used to condition dynamical systems on all manner of 
observed or desired qualitative dynamics, such as threshold behaviour, bifurcations, robustness, 
temporal ordering etc.. To rule out that a mathematical model can exhibit a certain dynamical 
behaviour will, however, require exhaustive numerical sampling of the parameter space; but 
coupled to ideas from probabilistic computing [27], our procedure lends itself to such investi- 
gations. Both for inference and design problems we foresee vast scope for applying this type of 
qualitative inference-based modelling. There is currently still a lack of understanding between 
the interplay of qualitative and quantitative features of dynamical systems J28| |; this becomes 
more pressing to address as the systems we are considering become more complicated and the 
data collected more detailed. Flexibility in parameter estimation — whether based on quali- 
tative or quantitative system features — will be an important feature for the analysis of such 
system, as well as the design of synthetic systems in engineering and synthetic biology. 

Methods 

Encoding dynamics through Lyapunov exponents 

Consider a continuous time dynamical system — similar results hold for the discrete case — 
described by, 

£ = (5) 

where / is an n-dimensional gradient field. To study the sensitivity of / to initial conditions, 
we consider the evolution of an initially orthonormal axes of n vectors, {ei, e 2 , e n }, in the 
tangent space at y . At time t, each e$ satisfies the linear equation, 

^ = Df(y t ) ■ ^ (6) 
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where Df(y t ) is the Jacobian of / evaluated along the orbit y t . Equations (2) and (3) describe 
the expansion/contraction of an n-dimensional ellipsoid in the tangent space at y t9 and we denote 
the average exponential rate of growth over all t of the ith principal axis of the ellipsoid as A$. 
The quantities, Ai > A 2 > ... > A n , are called the global Lyapunov exponents (LEs) of /. In 
particular, the sign of the maximal LE, Ai, determines the fate of almost all small perturbations 
to the system's state, and consequently, the nature of the underlying dynamical attractor. For 
Ai < 0, all small perturbations die out and trajectories that start sufficiently close to each 
other converge to the same stable fixed point in state-space; for Ai = 0, initially close orbits 
remain close but distinct, corresponding to oscillatory dynamics on a limit-cycle or torus (for 
tori, at least one other exponent must be zero); and finally for Ai > 0, small perturbations grow 
exponentially, and the system evolves chaotically within the folded space of a so-called "strange 
attractor" (for two or more positive definite LEs we speak of "hyperchaos"). 

In general, non-linear system equations and the asymptotic nature of the LEs precludes any 
analytic evaluation. Instead, various methods of numerical approximation of these quantities, 
both directly from ODE models and from time-series data p9] - [3T| have been developed. In 



this paper, Lyapunov spectra are calculated using a Python implementation of a method pro- 
posed by Benettin et.al. [32] and Shimada and Nagashima [33], (outlined in the Supplementary 
Information online) for inference of Lyapunov exponents when the differential equations are 
known. 

For each of the results presented in this article, we used LSODE to integrate the equa- 
tions initially for 1000 time steps, in order to overcome the transient dynamics. The Lyapunov 
exponents were then estimating over a further 10,000 points. The step size varied between 
0.01 for the Lorenz system and 0.5 for the Hesl model. The accuracy of our implementa- 
tion may be gauged from the sum of squared errors shown for the oscillation inference results 
where the maximal Lyapunov exponent should in theory be 0. For the Hesl system, limit 
cycles attractors were obtained with estimated | max(A^)| < 6 x 10" 3 , while for the electric 
ciricuit, the oscillations found had |max(A^)| < 3 x 10" 3 . Further, for the Lorenz system 
with typical parameter values (a, p, 0) = (10, 28, 8/3), we estimate the Lyapunov spectrum 
as (0.886, -4 x 10~ 3 , -15.2) as compared to the values (0.906, 0, -14.57) reported by J.C. 
Sprott j34). In this light, a conservative choice for 5 to i, the tolerance level above which we take 
an estimated maximal Lyapunov exponent to indicate the presence of chaos would be w 0.05 

Lyapunov spectrum driven parameter estimation 

Unlike in the case for linear systems, where identifying suitable parameters that produce ob- 
served or desired dynamics is trivial, inference for highly non-linear systems is far from straight- 
forward. Indeed, exact inferences are prohibitively expensive for even small systems, and so 



a host of different approximation methods have been proposed |20|[35}|36j. In our case, two 
further complications arise from using LEs to encode the desired behaviour. Firstly, the form 
of the mapping between model parameters and LEs is not closed, making methods that rely 
on an approximation of the estimation routine or its derivatives, such as the extended Kalman 
filter, difficult to apply. Secondly, LEs are significantly more expensive to compute than more 
traditional cost functions, ruling out the use of approaches such as particle filtering or sequen- 
tial Monte-Carlo methods that require extensive sampling of regions of parameter space and 
calculation of the corresponding LEs at each iteration. 



To overcome these challenges we exploit the efficiency and flexibility of the UKF [ 37-39 1, 



seeking here to infer the posterior distribution over parameters that give rise to the desired LEs. 
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Typically the UKF is applied for parameter estimation of a non-linear mapping g(-) from a 
sequence of noisy measurements, y k , of the true states, x k , at discrete times k = t l5 .., t N . A 
dynamical state-space model is defined, 

Ok = Ok-i + v k (7) 
Vk = 9{x k) e k )+u k (8) 

where u k -\ ~ N(Q, Q k ) represents the measurement noise, v k -i ~ N(0, R k ) is the artificial 
process noise driving the system, and <?(•) is the mapping for which parameters 9 k are to be 
inferred. The UKF (described in full below) is then characterised by the iterative application of 
a two step, predict and update, procedure. In the prediction step the current parameter estimate 
is perturbed by the driving process noise v k forming a priori estimates (which are conditional 
upon all but the current observation) for the parameter mean and covariance. These we denote 
as 6^ r and P^ r , respectively. The update step then updates the a priori statistics using the 
additional measurement, y k9 to form a posteriori estimates, 6 k ° and P k °. After all observations 
have been processed we arrive at the final parameter estimate, 6%° (with covariance Pf°). 

A crucial step in the algorithm is the propagation of the a priori parameter distribution 
statistics through the model, g(-). Assuming linearity of this transformation, a closed form 
optimal filter may be derived (known as the Kalman filter). However, this assumption would 
make the algorithm inappropriate for use with the highly non-linear systems and the choice 
of g(-) considered here. It is how the UKF copes with this challenge, namely its use of the 
"unscented transform", that makes it particularly suitable for our method of qualitative feature 
driven parameter estimation. 

The unscented transform is motivated by the idea that probability distributions are easier to 
approximate than highly non-linear functions J40| . In contrast to the Extended Kalman filter 
where non-linear state and transition functions are approximated by their linearised forms, the 
UKF defines a set, Q k , of "sigma-points" - deterministically sampled particles from the cur- 
rent posterior parameter distribution (given by Q v k _ x and P%° ), that along with corresponding 
weights, {c^ m ,^ c }&, completely capture its mean and covariance. The mean and covariance 
of the predicted observation, y k , may then be calculated to third order accuracy in the Taylor 
expansion, using the equations given below. Under the approximate assumption of Gaussian 
prior and posterior distributions (higher order moments may be captured if desired at the cost 
of computational efficiency), the deterministic and minimal sampling scheme at the heart of the 
filter requires relatively few LE evaluations at each iteration (2n p + 1, where n p is the number 
of parameters to be inferred). Further, the function that is the subject of the inference may be 



highly-nonlinear and can take any parametric form, such as a feed-forward neural network |4TJ, 
or as in our case, a routine for estimating the Lyapunov exponents of a model with a given 
parameter set. 

With [X]i denoting the i th column of the matrix X, the UKF algorithm for parameter esti- 
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mation is given by, 

Initialize: 

01° = E(6) 

P p ° = E(0 O - 6)00 - 6) T ) 
For each time point k — t 1? t N \ 
Prediction step: 

= Ci 

n r = p p k °-i+Qk-i 

Update step: 
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Various schemes for sigma-point selection exist including those for minimal set size, higher 
than third order accuracy and (as defined and used in this study) guaranteed positive-definiteness 
of the parameter covariance matrices |4^|42]-[44), which is necessary for the square roots ob- 
tained by Cholesky decomposition when calculating the sigma-points. The scaled sigma-point 
scheme thus proceeds as, 



[e fc ] = € 
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and parameters k, a and j3 may be chosen to control the positive definiteness of covariance 
matrices, spread of the sigma-points, and error in the kurtosis respectively. 

To apply to the UKF for qualitative inference, we amend the dynamical state-space model 

to, 



k 



Harget 



k -l + V k 

L(O k ,y ;f) +u k , 



(9) 
(10) 
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where L(-) maps parameters to the encoding of the dynamical behaviour (here a numerical rou- 
tine to calculate the Lyapunov spectrum), A targe t is a constant target vector of LEs, y denotes the 
initial conditions, and / is the dynamical system under investigation (with unknown parameter 
vector 0, considered as a hidden state of the system and not subject to temporal dynamics). To 
see how equations (4) and (5) fit the state- space model format for UKF parameter estimation it 
is helpful to consider the time series (A target , A targe t, A targe t, ...) as the "observed" data from which 
we learn the parameters of the non-linear mapping L(-). Our use of the UKF is characterised 
by a repeated comparison of the simulated dynamics for each sigma-point to the same (as spec- 
ified) desired dynamical behaviour. In this respect, we use the UKF as a smoother; there is 
no temporal ordering of the data supplied to the filter since all information about the observed 
(target) dynamics is given at each iteration. From an optimisation viewpoint, the filter aims to 
minimise the prediction-error function, 



E(0) = X>(0,3/o;/) - AtargetF^rMfK^O;/) " A target ], (11) 

1=1 

thus moving the parameters towards a set for which the system exhibits the desired dynamical 
regime. 

In the examples presented here, the measurement noise covariance is chosen as Q k = a/, 



with a G R and I the identity matrix, as suggested by E. Wan and R. Van Der Merwe [[451. 
We find that varying a over different orders of magnitude does not effect the ability to achieve 
qualitatively acceptable parameter regimes (results not shown); indeed it can be shown that 
(if the filter converges) a fixed diagonal measurement noise covariance matrix cancels out of 
the UKF parameter estimation algorithm [|45). However, we do find that the choice of a can 
influence the time taken to reach qualitatively acceptable parameter combinations, with a « 
0.01 performing well for the examples presented here. 

Non-diagonal entries of the process noise covariance, P k , are also fixed at zero, with each 
diagonal entry taken at the same order of magnitude as its corresponding initial model parame- 
ter choice. It is worth noting here that unlike in the state-estimation case, the "artificial" process 
noise, v k , has no physical interpretation. It may even be set to zero j46), though non-zero 
choices can help the algorithm skip out of non-optimal local minima (45j and allows the poste- 
rior to converge to non-point estimates, or converge at all in the case that the cross-covariance 
of the parameter prior and predicted data (Lyapunov exponents) remains non-zero. This can 
be seen by examining the UKF equations for converged P k °. Different methods exist for up- 
dating the process noise at each iteration of the filter to control the weight given to past and 
current observations J45] ], e.g. by annealing the covariance towards zero or by making use of 
the Robbins-Monro stochastic approximation scheme described in detail elsewhere |47). In the 
examples presented here, we find that keeping the value fixed gives the most reliable results, 
possibly reflecting the complex nature of the likelihood surface defined by our choice of g. 

Hes 1 Quantitative real-time PCR 

Dendritic cells (DC) were differentiated from bone marrow as described previously [48]. Rat 
Jgdl/humanFc fusion protein (R&D Systems) or human IgGl (Sigma Aldrich) (control sam- 
ples) were immobilised onto tissue culture plates (10 /xg/ml in PBS) overnight at 4C. DC were 
spun onto the plate and cells were harvested at the appropriate time. Total RNA was isolated 
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using the Absolutely RNA micro prep kit (Stratagene). cDNA was generated from 125 ng of 
total RNA using an archive kit (Applied Biosy stems). 1 fA of cDNA was used with PCR Mas- 
termix and TaqMan primer and probes (both Applied Biosystems) and analysed on an Applied 
Biosy stems 7500 PCR system. Cycle thresholds were normalised to 18S and calibrated to a 
PBS treated control sample for relative quantification. 



Computational Implementation 

All routines were implemented in Python using LSODE for integrating differential equations. 



ABC inference was performed using the ABC-SysBio package [ |49| |. Code is available from the 
authors upon request. 
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Appendix 



Lyapunov Exponents 

Lyapunov exponents describe the rate of separation of nearby trajectories and allow the pre- 
dictability of a system's future states to be quantified (see Figure 1 in the main paper). In our 
qualitative inference framework we exploit their ability to discriminate between qualitatively 
different orbit types, allowing us to drive the inference of parameters and initial conditions. 
Various algorithms exist for the estimation of these quantities both directly from ODE mod- 
els and from time-series data j29] - [3T| .For the results presented here, Lyapunov spectra were 



calculated using a Python implementation of a method proposed by Benettin et.al. p2| ] and 
Shimada and Nagashima (33), (outlined below) for inference of Lyapunov exponents when the 
differential equations are known. 



Estimating the Lyapunov spectrum of a differential equation model 

While analytic evaluations of Lyapunov spectra exist for certain special cases or simple sys- 
tems, generally applicable strategies must resort to numerical estimation techniques. A naive 
approach to studying a system's sensitivity to initial conditions would be to directly track the 
evolution, under /, of a set of initially close points and subsequently extract principal rates and 
directions of contraction/expansion. However, for the following reasons this turns out to be 
unfeasible: 

• The Lyapunov spectrum describes the average local rates of divergence of nearby trajec- 
tories. Under chaotic dynamics and given sufficient time, any deviation from an initial 
point, no matter how small, will grow too large to represent the local dynamics. 

• Finite errors in computer calculations and storage mean that every direction in state- space 
is contaminated by a component in the direction of the dominating max(A^). Hence 
any principal axis evolved through computer simulation will degenerate to align almost 
entirely along the direction of maximal expansion. 



A method employing Gram-Schmidt Re-orthonormalization (GSR) addresses these two is- 
sues [32, 33] (see Fig. SI). Here an initial TV-dimensional orthonormal axis, {ej, with 
origin at x , is chosen arbitrarily (e.g. the canonical basis) and integrated simultaneously to the 
initial condition. Whilst the trajectory from x is determined by /, each e\ is evolved according 
to Df(x t ). The first problem is thus avoided since the linearized equations approximate the real 
dynamics only at infinitesimally close points to x t . Linearity allows the application of Df(x t ) 
to finite magnitude vectors making up the principal axes. 

The second issue is resolved by a periodic application of GSR to re-orthonormalize the 
collapsing axes. Let {e*, e f N } be the set of vectors obtained by numerical integration of the 
principal axes at time t\ GSR defines a new orthonormal set {e\ , e^} as, 
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Note that the space spanned by the first k vectors is fixed under GSR, and so is free to "seek" 
(driven by numerical error) the /c-dimensional space with the highest rate of expansion. Thus, 
since the axes obtained by GSR are orthonormal, the fc-largest Lyapunov exponents may be 
obtained from the average rate of growth of the projection of the e\ onto e\. 




Figure 4: Lyapunov exponents. Lyapunov exponents characterise the long term evolution of the axes of an 
infinitesimal ball about an initial point in phase space. In the estimation method, an arbitrarily chosen orthonormal 
axes and the initial condition are evolved simultaneously via the linearised and true system equations respectively. 
GSR is periodically employed to correct for numerical corruption of the vector directions. Lyapunov exponents 
are calculated as the average rates of growth of the projections of the evolved axes vectors (blue) onto their re- 
orthonormalized versions (red). 



Constraining the inference 

Sometimes we may wish to constrain a search to particular regions of parameter space. In the 
context of modelling a real world process, this may be based upon the physical impossibility of 
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certain parameter combinations (e.g. negative chemical reaction rates). More generally we may 
wish to avoid 'badly behaving' regions of parameter space where the model is, for example, 
unbounded. Instead of constraining the unscented Kalman filter algorithm, we write a new 
observation function g* = g o p, where p maps the input parameters onto the region of interest. 
For example, in order to avoid negative chemical reaction rates, p may output the absolute value 
of the parameters (and the unchanged model). Note that the parameters inferred by the filter 
must then be interpreted in light of p. 

Model equations 

Below we state the governing equations for each of the systems considered. Details of the other 
models may be found in the references provided in the Results. For each system, a dot above a 
variable indicates the derivative with respect to time. 

Lorenz map 



X 


= <r(y 


— x) 


y 


= x(p 


-*) 


z 


= xy- 


-f3z, 



Chaotic electronic circuit 

x = y 

y = ay — x — z 
ez = b + y-c(e z -l), 

Hesl regulatory model 

M = -k deg M + l/(l + (P 2 /P ) h ) 
A = -k deg Pi + vM - fcxPx 
A = -k deg P 2 + fciPi, 

Four dimensional hyperchaotic system 

x\ = a(x2 — Xi) + x 2 x 3 

x 2 = b(xi + x 2 ) - xix 3 

X3 = —CX3 — 6X4 + x\x 2 

x\ = —dx 4 + fx 3 + X\x 2 
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